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1.0  INTRODUCTION 

Over  the  past  decade,  sufficient  advances  have  been  made  in  the  field  of  computational 
fluid  dynamics  (CFD)  to  provide  an  extremely  useful  tool  to  engineering  practitioners.  This 
tool  is  the  CFD  technology  which  currently  sees  extensive  use  in  a  large  variety  of  aerospace 
engineering  design  and  analysis.  The  present  report  deals  with  an  extension  of  this 
technology  which  is  immediately  applicable  to  the  study  of  flow  past  specific  components  of 
some  modern  aerospace  vehicles. 

Consider  the  wall/ramp  arrangement  shown  in  Fig.  1.  This  two-dimensional 
configuration  may  represent  a  simple  model  of  a  vehicle  control  surface.  A  ramp-induced 
shock-boundary-layer  interaction  results  when  this  model  encounters  a 
supersonic/hypersonic  airstream.  A  popular  computer  code  for  analyzing  this  type  of 
problem  is  documented  by  Hung  and  MacCormack  (Ref.  1)  and  by  MacCormack  (Ref.  2). 
This  code  (herein  referred  to  as  the  original  code)  solves  the  Navier-Stokes  equations  for 
two-dimensional  compressible  supersonic/hypersonic  laminar  flows  and  for  turbulent  flows 
which  are  adequately  described  by  an  algebraic  eddy  viscosity  turbulence  model. 

The  first  goal  of  the  present  effort  was  to  modify  the  original  code  to  apply  it  to  the 
cylinder-flare  configuration  illustrated  in  Fig.  2.  The  desire  to  compute 
supersonic/hypersonic  flow  over  cylinder-flare  configurations  is  motivated  by  the  following 
two  considerations.  First,  such  configurations  exist  in  real  aerospace  applications.  The 
second  and  perhaps  more  subtle  reason  is  related  to  the  difficulties  encountered  with 
comparing  results  from  the  original  two-dimensional  code  to  “two-dimensional” 
experimental  data.  The  “two-dimensional”  experimental  data  suffer  from  finite  span  effects 
which,  for  turbulent  flow,  may  be  important.  The  hollow  cylinder-flare  eliminates  these 
effects  from  the  experimental  data  so  that  less  uncertainty  exists  when  experimental  data  and 
the  appropriate  numerical  solution  are  compared.  Also,  it  remains  unclear  what  role,  if  any, 
three-dimensional  turbulence  plays  in  a  “two-dimensional”  problem.  Because  of  these 
considerations,  the  hollow  cylinder- flare  configuration  should  provide  a  better  shape  on 
which  to  base  the  design  of  new  turbulence  models. 

The  second  goal  of  the  present  effort  was  to  install  a  wall  mass  transfer  option  into  the 
computer  code  which  would  allow  treatment  of  the  hinge-line  boundary-layer  bleed  problem 
depicted  schematically  in  Fig.  3.  This  option  would  also  allow  treatment  of  the  boundary- 
layer  suction  problem  which  arises  frequently  in  aerospace  applications.  The  boundary-layer 
suction  technique  may  be  used,  for  example,  to  enhance  control  surface  effectiveness,  etc. 
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Figure  3.  Hinge-line  boundary-layer  bleed. 


Previous  investigators  have  directed  a  considerable  amount  of  effort  toward  the  study  of 
high-speed  shock-boundary-layer  interactions.  Of  primary  interest  for  the  present  study  are 
investigations  of  those  interactions  which  are  compression  corner-induced.  Laminar 
interactions  in  two  dimensions  have  been  studied  experimentally  by  Gray  and  Rhudy  (Ref.  3) 
and  theoretically  by  Bloy  and  Georgeff  (Ref.  4)  and  Carter  (Ref,  5),  for  example.  Turbulent 
interactions  in  two  dimensions  have  been  examined  experimentally  by  several  investigators 
including  Law  (Ref.  6),  Settles,  Vas,  and  Bogdonoff  (Ref.  7),  and  Gray  and  Rhudy  (Ref.  3). 
Investigators  of  theoretical  turbulent  interactions  include  Hung  and  MacCormack  (Ref.  1), 
Shang  and  Hankey  (Ref.  8),  and  Coakley,  Viegas,  and  Horstman  (Ref.  9).  Joint  theoretical- 
experimental  turbulent  flow  efforts  have  been  conducted  by  various  investigators  including 
Horstman,  Hung,  Settles,  Vas,  and  Bogdonoff  (Ref.  10).  Much  of  the  theoretical  effort  in 
turbulent  flows  has  been  devoted  to  the  testing  of  various  types  of  turbulence  models  in  an 
attempt  to  predict  with  accuracy  quantities  such  as  skin  friction  and  heat  transfer  for 
separated  flows. 

Literature  on  turbulent  interactions  for  axisymmetric  flows  is  not  as  abundant  as  for 
two-dimensional  flows.  Experimental  studies  are  documented  by  Kuehn  (Ref.  11),  Polak 
and  Kalivretenos  (Ref.  12),  Coleman  and  Stollery  (Ref.  13),  and  Roshko  and  Thomke  (Ref. 
14);  however,  numerical  studies  of  these  interactions  using  the  Navier-Stokes  equations  are 
scarce  in  the  literature.  It  is  typical  for  investigators  to  study  the  two-dimensional  problems 
first,  and,  once  these  are  understood,  to  extend  their  investigation  to  axisymmetric  and 
three-dimensional  problems.  Even  after  extensive  efforts,  two-dimensional  turbulent  shock- 
boundary-layer  interactions  for  supersonic/hypersonic  flows  are  not  yet  fully  understood. 
Therefore,  it  may  be  necessary  in  this  case  to  focus  first  on  the  axisymmetric  problem  due  to 
the  very  nature  of  turbulence.  At  least  part  of  the  difficulty  in  understanding  this  type  of 
flow  is  attributable  to  the  lack  of  understanding  of  the  magnitude  to  which  turbulence  levels 
are  influenced  by  three-dimensional  effects  such  as  the  finite  span  effects  previously 
mentioned.  This  is  not  to  say  that  a  concentrated  effort  on  the  axisymmetric  problem  will 
yield  all  the  answers,  but  rather  that  there  may  well  be  fewer  unknowns  with  which  to  deal. 
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The  additional  complication  of  wall  mass  transfer  for  flat  plates  has  received  extensive 
experimental  and  numerical  study;  however,  the  numerical  approaches  typically  involve  the 
boundary-layer  equations.  The  present  work  illustrates  how  the  wall  mass  transfer  idea  is 
applied  to  the  complete  Navier-Stokes  equations  for  compression  corner-induced  shock¬ 
boundary-layer  interactions. 

The  organization  of  the  remainder  of  this  report  is  now  described.  Section  2.0  deals  with 
the  governing  equations,  the  axisymmetric  assumption,  the  approximations  inherent  in  the 
turbulent  flow  equations,  and  an  assessment  of  these  approximations  for  some  of  the 
applications  of  interest.  This  section  also  contains  a  description  of  two  algebraic  eddy 
viscosity  turbulence  models  used  in  the  present  work.  Section  3.0  deals  with  the  numerical 
aspects  of  incorporating  the  axisymmetric  assumption  into  the  existing  computer  code. 
Specific  topics  addressed  include  the  finite-difference  grid,  the  time-splitting  approach,  and 
stability  of  the  resulting  source  term  operators.  A  description  of  the  analysis  and  numerical 
application  involved  in  generalizing  the  boundary  conditions  to  include  the  effectss  of  wall 
mass  transfer  is  included  in  Section  4.0.  Section  5.0  contains  comparisons  of  numerical 
results  with  experimental  data  which  verify  the  resulting  computer  code  for  a  certain  class  of 
problems  and  illustrate  its  deficiencies  for  another  class  of  problems.  A  discussion  of  the 
code's  potential,  as  well  as  its  deficiencies,  is  also  presented.  Section  6.0  contains  concluding 
remarks. 


2.0  GOVERNING  EQUATIONS 
2.1  THREE-DIMENSIONAL  NAVIER-STOKES  EQUATIONS 


The  fully  three-dimensional  time-dependent  Navier-Stokes  equations  are  given  in 
conservation-law  form  by 


dd>  d&  df  dg 
— —  +  —  +  —  +  -5 
tit  <9x  dy  dz 


=  0 


(1) 


where 


e  =■  e  ■  —  e  -  7  p  —  & 

i  v  *  i  Lv  ’  g  -  gj  gv 

and  the  barred  quantities  indicate  an  association  with  the  Cartesian  reference  frame.  The 
subscripts  i  and  v  denote  inviscid  and  viscous  contributions,  respectively.  The  quantities  e,, 
fi,  and  a  are  given  by 
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The  solution  vector  0  is  given  by 


In  these  equations,  /t  is  the  dynamic  or  first  coefficient  of  viscosity,  X  =  x  -2/3  n  is  the 
second  coefficient  of  viscosity*  (actually,  as  is  typically  done  in  practice,  x,  the  bulk 
viscosity  coefficient,*  is  taken  to  be  zero),  k  is  the  thermal  conductivity  of  the  fluid,  e  is  the 
total  energy  per  unit  volume  of  the  fluid  and  is  related  to  pressure,  p,  density,  q,  and  (x,  y,  z) 
velocity  components,  (u,  v,  w),  by  the  equation 

*  =  +  2) 
y- 1  2 

assuming  the  fluid  behaves  according  to  the  perfect  gas  law  p  =  (7-l)e  e"  where  e  is  the 
specific  internal  energy  of  the  fluid.  In  addition,  T  is  the  static  temperature  of  the  fluid,  q 

A  A  A 

is  the  velocity  vector,  (i,  j,  k)  are  the  Cartesian  unit  vectors  in  the  (x,  y,  z)  directions, 
respectively,  and  V  is  the  gradient  operator. 


*The  rotation  and  definitions  advocated  by  Rosenhead  (Ref.  15)  and  used  by  Owczarek  (Ref.  16) 
are  employed  here.  These  differ  from  those  used  by  Lamb  (Ref.  17)  and  Yuan  (Ref.  18), 
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2.2  THE  AXISYMMETRIC  ASSUMPTION 

To  treat  the  axisymmetric  flow  case  it  is  necessary  to  introduce  a  more  general  set  of 
independent  variables.  Let 

r  -  t 

f  =  * 

r]  =  7?(y,  z) 

£  - 


Equation  (1)  may  be  rewritten  in  terms  of  these  independent  variables  to  yield 


d_ 

ec 


-  [(’y)  -(^)JT-  [(’.)/  (iz)Ji-  ■  0 


(2) 


These  independent  variables  are  made  to  look  like  cylindrical  coordinates  by  introducing  the 
rotation 

y  =  tjcos£ 


7.  =  )jsin£ 

from  which  the  following  relations  are  obtained: 

sin  £ 


J 


r)y  =  cos£  .  r/z  =  sin£  ,  £y  = 


n 

cos  £ 


("y)„  -  °  •  (<y)(  ’  ~i 


i  * 


(’0,  ■ 0  ■  (<-)4  - 


sin  4 


- -  =  C 


The  use  of  these  relations  in  Eq.  (2)  produces  a  repeated  occurrence  of  the  groupings  {ijyf  + 
ijzg)  and  (?jyg  -  Tjzf).  Since  the  cylindrical  (£,  tj,  f )  velocity  components  (u,  v,  w)  are  related  to 
the  Cartesian  components  by 

U  =  U 

v  =  r)yV  tjzw 
w  =  7?y-vv-7jz7 
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it  is  a  simple  geometric  exercise  to  show  that  replacing  the  third  component  of  Eq.  (2)  by  ijz 
times  the  fourth  component  plus  r)y  times  the  third  component,  and  replacing  the  fourth 
component  of  Eq.  (2)  by  times  the  fourth  component  minus  jj,  times  the  third  component 
yields  the  Navier-Stokes  equations  in  cylindrical  coordinates: 


where 


dd>  dc  d  f  d  £  £  . 

dr  dr)  d£ 


A 

<f>  =  (p  ,  pu  ,  pv  ,  pw  ,  e) 


A 

e 


A 

ei 


A 

■  c. 


A 

f 


AAA 
8  “  §  i  8  v 


A 

f: 


A 

f. . 


are  given  by 


"The  superscript  T  denotes  the  transpose  operation. 
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where  1$,  i,,  if  are  the  cylindrical  unit  vectors  in  the  direction  indicated  by  the  subscripts. 
These  equations  look  like  their  counterparts  in  the  Cartesian  reference  frame  except  for  the 

A  a 

appearance  of  the  source  terms  hi  and  hv  and  the  terms  encircled  by  dotted  lines  in  the 

^  A 

expressions  for  fv  and  gv. 


If  the  problems  of  interest  are  restricted  to  nonspinning  bodies  of  revolution  at  zero 
incidence,  then  the  meridional  velocity  and  its  gradient  are  zero.  In  this  case  the 
{■-momentum  equation  is  trivially  satisfied  and  the  remaining  system  of  equations  is  given  by 


d6  d*_ 
<5«  dx 


6M[_ 

<9y 


Oi  =  0 


(3) 


where  the  independent  variables  (t,  x,  y)  have  been  re-introduced  for  readability  with  the 
obvious  associations 

T  1  >  £  x  •  v  y  » 


The  quantities  <f>,  e,  f,  and  h  are  given  by  <f>  =  (c,  eu,  qv,  e)T;  e  =  ej-  ev;  f  =  f;  -  fvt  and  h  = 
h;  -  hv  where 
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e 


f 


V 


h 


v 


1 


y 


(4) 


(5) 


(6) 


These  equations  appear  to  be  identical  to  their  two-dimensional  Cartesian  counterparts 
except  for  the  source  term,  h.  There  is  also  the  hidden  difference  in  the  velocity  divergence 
term,  V.q,  which  for  the  axisymmetric  case,  is  given  by 


V  ■  q  = 


<3  u 
dx 


(9v  V 


<3y  y 


(7) 


As  is  demonstrated  in  Section  3.0,  this  similarity  between  the  two-dimensional  Cartesian  and 
the  axisymmetric  equations  allows  the  axisymmetric  assumption  to  be  incorporated  into  the 
original  two-dimensional  computer  code  as  a  pair  of  source  term  operators  in  the  already 
present  time-split  operator  sequence. 
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2.3  TURBULENCE  MODELING 
2.3.1  Concepts  and  Assumptions 

It  is  generally  accepted  that  the  Navier-Stokes  equations  are  applicable  to  both  laminar 
and  turbulent  flows  (Ref.  19).  It  is  unfortunate,  however,  that  the  small  significant  length 
scales  present  in  turbulent  flows  render  the  numerical  simulation  of  such  flows  impractical 
with  regard  to  resolving  these  scales  on  current  computer  architecture.  It  is  therefore 
necessary  to  employ  some  type  of  modeling  procedure  to  compute  turbulent  flows. 

The  mass-weighted  averaging  technique  usually  associated  with  compressible  flows  is 
used  in  the  present  approach.  A  description  of  this  technique  is  given  by  Marvin  (Ref.  19). 
The  key  feature  of  this  method  is  that  there  are  fewer  terms  representing  the  so-called 
turbulent  Reynolds  stresses  and  heat  fluxes  than  for  the  Reynolds  averaging  method.  The 
modeling  problem  may  thus  be  reduced  to  a  problem  of  representing  these  stresses  and  heat 
fluxes  in  terms  of  the  mean  and  mass-weighted  flow  variables.  One  approach  is  to  use  the 
classical  Boussinesq  assumption,  which  provides  a  simple  means  to  simulate  the  effect  of 
turbulence  through  an  eddy  viscosity  formulation  where  the  n  appearing  in  Eqs.  (4)  through 
(6)  is  simply  taken  to  be  the  sum  of  a  laminar  contribution,  m,  and  a  turbulent  contribution, 
Ht.  The  laminar  viscosity  is  assumed  to  be  governed  by  Sutherland’s  Law  given  by 


n  =  2.27(10) 


-s 


...  3/  2 


]’  +  198.6 


where  the  static  temperature,  T,  is  in  degrees  Rankine.  The  turbulent  viscosity  is  determined 
from  a  two-layer  model  to  be  discussed  later.  A  treatment  of  the  heat  flux  terms  of  Eqs.  (4) 
through  (6)  uses  the  definition  of  Prandtl  number,  Prf,  and  the  introduction  of  a  turbulent 
Prandtl  number,  Prt 


Pr<  = 


Pr. 


where  cp  is  the  specific  heat  at  constant  pressure.  The  conductivity,  k,  in  Eqs.  (4)  through  (6) 
is  then  given  by 

k-k, 

"V-i  »>) 

where  both  Pr?  and  Prt  are  assumed  to  have  constant  values  of  0.72  and  0.9,  respectively. 


As  mentioned,  the  turbulent  viscosity  is  modeled  with  a  two-layer  concept.  The  basic 
inner  layer  model  uses  Prandtl’s  mixing-length  hypothesis  and  is  given  by 
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(fit)  =  pt2  (y)G 

v  'inner 


(8) 


where  g  is  the  local  mean  density,  0(y)  =  k  y  D  where  k  is  a  constant  usually  taken  to  be  0.4, 
and  D  is  the  van  Driest  damping  factor  given  by 


D  = 


1  - 


e 


+ 


f 


where  yx  =  y  ^ Pt  and  A+  may  be  used  to  incorporate  the  effects  of  streamwise 

pressure  gradient  and  wall  mass  transfer  similar  to  its  use  in  incompressible  flows  as 
indicated  by  Marvin  (Ref.  19).  The  subscript  r  appearing  in  the  definition  of  y+  denotes  a 
reference  condition  evaluated  at  the  wall  point.  G  denotes  some  measure  of  the  local  mass- 
weighted  velocity  gradient. 


The  outer  layer  model  is  given  by 

(f*i)  =  CypQ(y)  (9) 

'  f  outer 

where  C  is  a  constant,  y  is  a  reference  distance  normal  to  the  wall,  and  Q(y)  is  a  Klebanoff- 
type  intermittency  factor  given  by 

Q(y>  =  1  +  5.5  (by)6]  _1  (10) 


where  b  is  a  constant. 

in  practice,  fa)inner  is  used  out  to  the  first  point  where  it  exceeds  (/Pouter*  and  On)oui<r  is 
used  from  this  point  out  to  the  edge  of  the  turbulent  layer. 

At  this  point  a  brief  summary  is  given  of  the  basic  turbulence-related  assumptions 
inherent  in  the  mass-weighted  averaged  Navier-Stokes  equations  used  in  the  present  work: 

1 .  It  is  assumed  that  the  goal  is  only  to  determine  the  mass-weighted  and/or  mean 
variables,  since  the  averaging  process  precludes  the  reconstruction  of  the  original 
dependent  variables  (i.e.,  a  loss  of  information  is  accepted). 

2.  Reynolds  stresses  and  turbulent  heat  flux  have  the  same  form  as  the 
conventional  laminar  stress  tensor  and  laminar  heat  flux,  respectively. 
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3.  The  laminar  and  turbulent  Prandtl  numbers  are  assumed  to  be  constant  and 
their  values  are  taken  to  be  0.72  and  0.9,  respectively. 

4.  The  eddy  viscosity  may  be  represented  by  a  simple  two-layer  algebraic  model  of 
the  type  previously  described. 

Although  all  of  these  assumptions  limit  the  applicability  of  the  given  equations,  probably 
the  most  restrictive  of  all  are  assumptions  2  and  4.  Efforts  are  underway  by  various 
investigators  to  relax  these  two  assumptions  somewhat  by  developing  more  rigorous  and 
more  complex  methods  to  evaluate  the  Reynolds  stresses  and  turbulent  heat  flux  terms. 
These  methods  attempt  to  model  more  fundamental  properties  of  the  turbulent  flow  and  are 
referred  to  as  multi-equation  models. 

2.3.2  Current  Model 

The  model  existing  within  the  original  computer  code  is  referred  to  as  the  current  model. 
This  model  is  styled  after  one  by  Cebeci,  Smith,  and  Mosinskis  (Ref.  20),  with  G  in  Eq.  (8) 
given  by 


dv 


where  u  is  the  component  of  mass-weighted  velocity  normal  to  the  y  direction  and  y  is 
measured  normal  to  the  wall.  C,  y,  and  b  in  Eqs.  (9)  and  (10)  are  given  by 

C  =  Kul} 

b  =  I 

S 

where  &  is  the  boundary-layer  thickness,  ue  is  the  velocity  at  the  boundary-layer  edge,  and  K 
is  an  empirical  constant  set  at  0.0168.  Unfortunately,  this  model  requires  that  a  boundary- 
layer  edge  be  clearly  distinguishable,  which  is  not  the  case  for  many  flows  with  shock¬ 
boundary-  layer  interactions,  including  separated  flows. 

2.3.3  Alternate  Model 


An  alternate  model  is  one  developed  by  Baldwin  and  Lomax  (Ref.  21),  with  G  in  Eq.  (8) 
given  by 
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and  C,  y,  and  b  in  Eqs.  (9)  and  (10)  given  by 

C  =  KCcp 

y  =  ^  WAKE 
b  =  CKLEB/yMAX 

Fwake  is  the  minimum  of  yMAX  Fmax  and  CwxyMAX  U2dif/Fmax  where  FMAx  is  the 
maximum  value  of  the  function 

F(y)  =  y|aj|D 

and  yMAX  is  the  location  at  which  FMax  occurs.  Udif  is  defined  by 

LDIF  =  (>/u2  +  V2  )MAX  -  (V“2  +  V2  )mIN 

and  the  constants,  CCp,  CKLeb,  and  Cwk  are  set  by  Baldwin  and  Lomax  to  the  values 

Ccp  =  1.6 

CKLEB  =  °-3 
CWK  =  0.25 

Note  that  this  model  does  not  require  the  determination  of  a  boundary-layer  edge. 

Many  variations  exist  in  the  two  models  just  described  but  the  worthiness  of  any  of  these 
models  is  a  function  of  the  particular  problem  being  solved. 

3.0  NUMERICAL  INCORPORATION  OF  THE  AXISYMMETRIC  ASSUMPTION 

3.1  FINITE-DIFFERENCE  GRID  AND  BOUNDARY  CONDITIONS 

As  mentioned  in  the  introduction  and  illustrated  in  Fig.  1 ,  the  original  computer  code 
was  designed  for  a  specific  class  of  problems  —  namely  two-dimensional  wall-ramp 
configurations.  The  purpose  of  this  section  is  to  illustate  how  an  extension  was  made  to 
allow  treatment  of  the  cylinder-flare  configuration  shown  in  Fig.  2.  The  basic  physical 
domain  appears  identical  to  its  two-dimensional  counterpart  since  only  one  plane  of  the 
axisymmetric  problem  is  required.  The  finite-difference  grid  depicted  in  Fig.  4,  therefore,  is 
identical  to  that  documented  by  Hung  and  MacCormack  (Ref.  1).  This  is  a  double  layer  grid 
with  the  first  layer  indicated  by  the  height,  HF.  The  grid  lines  in  this  layer  are  highly 
clustered  near  the  wall.  This  layer  is  called  the  fine  mesh  and  is  intended  to  contain  all  of  the 
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Figure  4.  Representative  finite-difference  grid 
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viscous  interaction  of  the  problem  including  the  boundary  layer  and  the  separation  zone,  if 
one  exists.  The  outer  layer,  or  coarse  mesh,  extends  from  the  outer  edge  of  the  fine  mesh  to  a 
height,  H,  above  the  wall.  This  layer  is  intended  to  contain  the  inviscid  features  of  the  flow 
such  as  leading  edge  shocks  or  bow  shocks.  The  outer  edge  of  this  layer  is  presumed  to  be  in 
the  free-stream  flow  region.  The  grid  point  indices  are  (i,  j)  as  indicated  in  the  figure. 


One  final  point  to  be  made  with  reference  to  Fig.  4  concerns  the  application  of  boundary 
conditions  on  this  grid.  The  upper  boundary  conditions  are  fixed  at  the  free-stream  flow 
values.  The  left  boundary  condition  is  fixed  at  some  known  inflow  solution.  The  right 
boundary  requires  some  type  of  simulation  of  an  outflow  condition.  This  is  accomplished  by 
a  simple  extrapolation  of  the  solution  along  j  =  constant  grid  lines.  The  lower  grid 
boundary  is  used  as  a  reflection  line  to  simulate  the  existence  of  a  solid  boundary.  The 
problem  here  is  to  determine  the  primitive  variables  (p,  T,  u,  v)  at  points  (i,  1)  based  on 
known  or  assumed  properties  of  viscous  flow  at  solid  boundaries.  With  reference  again  to 
Fig.  4,  points  ahead  of  the  solid  wall  (if  there  are  any)  where  i  <  if.e.  are  treated  very  simply 
by 

Pi.  1  “  Pi, 2 


TU  =  Ti,2 


ui.l  -  “i.2 


(ID 


The  points  where  i  2s  i*.e.  are  handled  by  applying  the  no-slip  velocity  condition 


uU  =  -Ui.2 

vi,l  =  -Vi,2 

(12) 

and  by  assuming  either  a  fixed  wall  temperature,  Tw,  such  that 

Ti.l  =  2TW-Ti,2 

or  an  adiabatic  wall  condition  such  that 

(13) 

Tu  -  [l— WT‘  — I.* 

(14) 

where  f(4> j)  =  sin  <j>-t  cos<fr  (yi.2-yi,iV(xi*xi-i)  aild  4>i  is  given  by 

(  0  for  i  <  ic 

(15) 

9i  ~  { 

(  8  for  i  >  ic 
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The  remaining  variable  pj,]  is  determined  by  solving  a  finite-difference  approximation  to 
the  steady  form  of  the  y-momentum  equation*  applied  at  the  solid  boundary.  This  equation 
is  obtained  by  performing  a  coordinate  mapping  on  Eq.  (3)  so  the  new  computational 
coordinate  lines  are  as  indicated  in  Fig.  4.  This  is  accomplished  by  introducing  the 
computational  coordinates  (r,  tj)  defined  as 


T  -  t 


f  =  * 


>'  -  >  D 


for  x  <  : 


1  = 


y  —  —  (x— xr)tan0  for  x  >.5^, 

Application  of  this  mapping  to  Eq.  (3)  yields 


where 


where 


d<f>  do 


:)  -  h 


0 


i  0  for  x  <  xc 

=  )  or  =  — tan0(x) 

i  —tan  9  for  s  >  xf 


V  0  for  x  <  \r 

0(x)  = 

jd  for  x  >  xr 


(16) 


The  pressure  is  then  obtained  by  setting  0T  =  0  and  solving  the  third  component  of  Eq.  (16) 
for  p,.  The  result  is 


=  [-pv(v^xul-,(v  +  ,xU)7-Mxvf+,(l  ^)vv  +  Aa£- A(v-^xu)^] 

77  i  ’  i  ^  07) 

+  | ’-PBV  +  pi?xv7J  -  ^(Av)^  +  —  pv 2  J  —  —  [(A  -  2p)vJ 

This  equation  has  the  form  p,  =  right-hand  side  (RHS)  where  RHS  is  evaluated  at  the  solid 
boundary.  If  p,  is  approximated  by  a  central  difference,  the  quantity  p,,,  is  obtained  from 


*The  use  of  the  steady  form  of  the  normal  momentum  equation  is  valid  for  unsteady 
flows  since  the  d/dt  term  involves  only  the  time  variation  of  the  wall-normal  mass  flux.  This 
variation  is  zero  for  time-variant  wall  mass  transfer. 


22 


AEDC-TR-80-24 


Pi,i  =  Pi.2-(y2-yi)RHS  (18) 

It  should  be  noted  that  for  the  specific  case  where  u  =  v  =  0  at  the  solid  boundary  (no-slip 
and  no  mass  transfer),  many  of  the  terms  in  RHS  vanish. 

3.2  TIME-SPLITTING  APPROACH 

The  time-splitting  approach  of  MacCormack  (Ref.  2),  also  known  as  the  method  of 
fractional  steps,  is  used  in  the  original  computer  code  to  solve  the  two-dimensional 
counterpart  of  Eq.  (16). 

A  complex  multidimensional  partial  differential  equation  may  be  split  into  a  sequence  of 
simpler  one-dimensional  equations,  each  of  which  is  then  approximated  by  finite- 
differences.  Then  a  simple  sequence  of  various  one-dimensional  operators  yields  successive 
candidate  solutions  of  which  the  final  one  approximates  the  solution  to  the 
multidimensional  equation.  This  allows  consideration  of  various  terms  or  groups  of  terms 
separately.  For  example,  consider  the  two-dimensional  equation  in  computational 
coordinates 

d$  3E  dF 

0 

where  $  =  0,  E  =  e,  F  =  f  +  jjxe  and  0,  e,  f  are  the  same  as  those  appearing  in  Eq.  (3).  The 
solution  to  this  equation  may  be  obtained  by  applying  the  following  sequence  of  operators: 

^  e2 

Oa  =  +0n-  —  E*  ) 

2  \  2  5/ 


0b  =  $>a-ArFa 

=  I  ^>b  +  S>a  -  \T  fJ  ) 

$n+  1  =  <J>b  _  —  Eb 
2  C 

ijjn+1  =  1  +  1  +  ob  -  —  E1+  1  'l 

2  V  2  f  / 
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so  that  in  operational  notation 

®”+l  ■  (TO 

The  operator  Lf  is  applied  twice  to  make  the  sequence  of  operators  symmetric  to  provide 
formal  second  order  accuracy  in  t.  It  is  equally  valid  to  apply  L,  twice  and  Lf  once,  but  it  is 
more  efficient  to  apply  least  often  the  operator  requiring  the  most  computational  effort. 
Unfortunately,  this  decision  may  not  be  optimal  if  stability  considerations  dictate  that  the 
operator  with  the  more  severely  restricted  stability  condition  be  applied  over  the  smallest 
step  size  interval  (At/2).  This  is  precisely  the  sequence  used  for  the  coarse  mesh  in  the 
original  computer  code.  For  the  fine  mesh,  however,  the  L„  operator  is  broken  up  into  a 
hyperbolic  part  involving  only  inviscid  terms  and  a  parabolic  part  involving  the  viscous 
contribution.  The  resulting  sequence  is  given  by 

*■*'  -  Li(f)L$(f)Lf(M  L$  (fK  (t)*’  (20) 


The  operator  Ljj  uses  an  approximation  to  the  method  of  characteristics  to  solve  the  equation 


d~  dy 


0 


(21) 


where  Fi  =  fj  +  tjx  e,.  It  is  shown  in  Section  4.0  that  this  approximation  is  troublesome  when 
wall  mass  transfer  is  considered.  The  operator  L£  solves  the  equation 


<M>  <3FV 

dr  <?JJ 


(F>  -  <, *  *,«,) 


with  an  implicit  technique  which  takes  advantage  of  the  parabolic  nature  of  the  equation. 


Of  course,  the  advantage  of  using  the  time-splitting  technique  is  that  a  less  severe 
stability  limitation  on  the  stepsize,  At,  is  realized  and  the  stability  analyses  are  simplified 
because  all  operators  are  one-dimensional. 


Since  the  only  difference  between  the  equations  just  described  and  Eq.  (16)  is  the  source 
term,  h,  and  the  extra  source  in  the  V.q  terms  of  ev  and  fv  [see  Eq.  (7)1,  a  simple  method  of 
incorporating  these  terms  into  the  existing  code  is  to  add  another  operator  to  the  time-split 
sequence  just  described.  This  operator  solves  the  equation 

-  =  0  (22) 

dr 
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where  5  is  the  contribution  due  to  the  source  in  V.q.  Let  this  operator  be  denoted  by  Lh6  so 
that  4»m  + 1  =  Lhi  where  the  superscript  m  indicates  a  fractional  step  operation.  Equation 
(22)  could  be  solved  in  a  number  of  different  ways;  however,  it  is  most  convenient  and 
efficient  to  introduce  the  primitive  variable  vector  as  the  new  dependent  variable  set.  The 
equation  is  then  given  by 

—  4-  =  0  (23) 

dr 

where  ^  =  (e,  u,  v,  p)T  and  h$,  h”  are  given  by 


ix  V  -  7 IAv), 


2\v  (VjJ  +  7,x  u^)  t  kT^ 


and  v  is  the  kinematic  viscosity,  h/q.  The  relation  between  <h  and  may  be  expressed  as 

$  =  (26) 

where  M  is  a  simple  nonlinear  algebraic  operator.  In  addition,  the  solution  to  Eq.  (23)  may 
be  expressed  as  ,  , 

where  U  is  an  appropriate  finite-difference  operator.  Since  a  clear  distinction  is  apparent 
between  terms  involving  d/drj  operations,  d/d£  operations,  and  algebraic  operations  in  Eqs. 
(24)  and  (25),  the  operator  is  split  into  two  parts  so  that 

=  L| l* 
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where  the  d/dy  terms  are  contained  in  L>  and  the  rest  of  the  terms  are  contained  in  L*.  The 
operator  Ljf  solves  the  equation 


=  0 


(27) 


and  solves 
v 


Sip 

57 


+  1.7  =  0 


(28) 


The  fractional  step  from  mtom+1  is  thus  expressed  as 


It 


III  \ 


=  i.|  lJ 


which  is  written  with  the  use  of  Eq.  (26)  as 


Premultiplying  by  M  yields 


W'J  (tm  +  1  =,  1,^  M'1 


<t>:" +  1  =  vii  ^  n_J  om 


Thus,  for  the  fine  grid,  the  entire  operator  sequence  for  one  symmetric  step  is  given  with 
reference  to  Eq.  (20)  as 

(29) 

An  extension  of  Eq.  (19)  for  the  coarse  grid  is  given  by 

'  LKr^7)s^"-'^?(r),-/(r)“-MrK  ,30) 

The  numerical  details  of  the  operators  Ljf  and  L*  are  now  presented. 

A  closer  look  at  Eq.  (27)  reveals  a  useful  structure  of  this  system  which  is  exploited  in  the 
development  of  the  operator  Ljf.  Substitution  of  Eq.  (24)  into  Eq.  (27)  yields  a  partial 
uncoupling  of  the  system  of  equations.  In  particular,  the  first  and  third  components  of  this 
system  are  written  as 

!  (31) 

(2/t  +  A)  v 
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If  it  is  assumed  that  p  and  X  are  time  invariant,  these  equations  represent  a  system  of  first 
order  ordinary  differential  equations  in  the  variables  q  and  v.  The  second  component  of  Eq. 
(27)  is  given  by 

-  »[p  <32> 

This  equation  is  also  an  ordinary  differential  equation  and  is  solved  for  u  very  easily  since 
the  entire  RHS  is  known  once  the  system  represented  by  Eq.  (31)  is  solved  for  q  and  v.  The 
fourth  component  of  Eq.  (27)  is 


Again,  the  RHS  is  known  and  the  coefficient  of  p  is  a  known  function  of  r  (known 
numerically),  so  another  ordinary  differential  equation  results  —  this  time  for  the  pressure, 
p.  Equations  (31)  through  (33)  may  be  easily  solved  numerically  and  variations  in  p  and  X  are 
accounted  for  by  using  a  predictor-corrector  algorithm.  That  is,  the  equations  are  solved  as 
described  with  current  p  and  X  values.  The  p  and  X  values  are  then  updated  based  on  this 
solution  and  the  corrector  step  is  taken.  Central  differences  are  used  for  all  spatial 
derivatives.  This  is  a  rapid  second-order  accurate  procedure  which  describes  the  operator 


The  operator  Ljj'  which  solves  Eq.  (28)  is  now  described.  The  first  component  of  this 
vector  equation  is 


which  is  trivially  solved  at  each  point  with  q  =  constant  during  this  step.  The  remaining 
three  equations  are  solved  using  MacCormack’s  standard  one-dimensional  scheme.  These 
three  equations  are 


+  -  o 

— - -  f  2fiv  +-  (Xv)  1.0 

dr  py  L  71  7  J 


Let  the  superscript  p  denote  the  predicted  value  and  the  superscript  c  denote  the  corrected 
value,  with  no  superscript  indicating  the  current  value.  The  numerical  scheme  for  solving 
these  equations  is  then  given  as 
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This  scheme  identifies  the  operator  Ljj\ 

At  this  point  the  two  additional  operators  and  A  have  been  described;  however,  their 
effect  upon  the  stability  of  the  entire  sequence  of  operators  must  be  determined. 

3.3  STABILITY 


A  sufficient  condition  for  the  stability  of  the  entire  sequence  of  operators  is  that  each 
operator  in  the  sequence  be  stable.  Therefore,  it  is  enough  to  determine  the  conditions 
required  for  each  of  the  operators  Lj?  and  to  be  stable  since  MacCormack  (Ref.  2)  has 
given  stability  conditions  for  the  other  operators  in  the  sequence. 


For  the  operator  L^,  consider  a  slightly  different  form  of  Eq.  (27)  written  as 


di/r  ^  Ir 

+  '  d( 


+  s 


where  A  is  a  Jacobian  matrix  given  by 


o 


0 


A  = 


0 


0  0  0  \ 

o  -  —{v  +  v')  o  ; 

Y 

0  0  0 


1  2(y— l)Av 

\  0  *~y -  0 


o 


/ 


with  v '  =  \/q  and  s  is  a  collection  of  nonderivative  terms*.  The  eigenvalues  of  the  matrix  A 
indicate  the  character  of  the  equations.  These  eigenvalues,  denoted  by  a,  are  all  zero.  Since 
the  well-known  CFL  (Courant,  Friedrichs,  Lewy)  stability  condition  states  that  the 
maximum  allowable  step  size.  At,  depends  on  the  maximum  eigenvalue  of  the  Jacobian 
matrix  through  the  relation 


this  clearly  implies  that  any  At,  which  satisfies  the  stability  requirement  of  the  other 
operators  will  satisfy  the  requirement  for  the  operation.  This  result  is  consistent  with  the 
fact  that  this  system  reduces  to  a  group  of  ordinary  differential  equations. 


♦Spatial  variations  in  the  second  coefficient  of  viscosity,  X,  are  neglected  in  this  analysis. 
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Now  consider  the  following  form  of  Eq.  (28)  for  the  Ljj'  operator: 


rX b  .  d<Jj 

—  +  A  — 

dr  f?7) 


H  S 


0 


where  A  is  given  by 


where  a2  =  Yp/7.  The  eigenvalues  of  this  matrix  are 


V  1  ^ 

CT|  =  0  ,  On  — - ,  ffn  - - (2l/  +  1/'  )  ,  n ,  —  -  - 

1  c  y  0  y  *■  ypr- 

In  practice,  X  =  -2/3  ji  so  that  v'  =  -2/3  v  and  cr3  =  -4/3  v/y.  Therefore  |<r3 1  >  \o2  | .  Also, 
k/gcv  =  yv/Pr  where  Pr  is  the  Prandtl  number.  Thus,  04  becomes  —  7p/(yPr).  The 
maximum  magnitude,  |a|ma!t,  is  then  given  by 


*;[“(;  £-H)] 
-(£  •  0"(?) 


The  quantity  7/Pr  will  typically  be  Larger  than  4/3  so  the  maximum  At  due  only  to  this 
operator  is  restricted  by  the  inequality 


I,! 


It  has  been  found  in  practice  that  this  stability  condition  is  typically  several  orders  of 
magnitude  less  restrictive  than  the  condition  which  already  exists,  due  to  the  other  operators 
in  the  sequence.  Therefore  the  present  method  of  introducing  the  axisymmetric  source  terms 
into  the  existing  computer  code  is  a  very  effective  one. 
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In  addition  to  generalizing  the  computer  code  for  application  to  axisymmetric  flow 
problems,  a  portion  of  the  current  effort  involves  modification  of  the  boundary  conditions 
to  allow  simulation  of  wall  mass  transfer.  This  the  subject  of  the  next  section. 

4.0  WALL  MASS  TRANSFER 
4.1  WALL  BOUNDARY  CONDITIONS 

Some  motivation  for  desiring  an  option  to  include  effects  of  wall  mass  transfer  is  given  in 
Section  1.0.  This  option  normally  requires  a  simple  modification  to  the  conventional 
conditions  applied  at  a  solid  boundary.  Unfortunately,  an  additional  modification  is 
required  in  the  present  application  because  of  an  implicit  assumption  in  the  original 
development  of  the  code  regarding  the  magnitude  of  the  wall-normal  flow  velocity  in  the 
fine  mesh.  This  modification  is  considered  a  major  one  and  will  be  discussed  following  the 
description  of  the  new  wall  boundary  condition  procedures. 

Let  mj  be  the  normal  mass  flux  to  be  specified  at  the  ith  point  along  the  solid  boundary 
grid  line  (see  Fig.  4).  Then 

mi  “  Pi,wVi.w 

where  Vi7W  is  the  flow  velocity  normal  to  this  boundary  and  q\^  indicates  the  wall  density  at 
this  point.  Then  the  actual  boundary  condition  procedure  is  similar  to  that  described  in 
Section  3.0  by  Eqs.  (1 1)  through  (15).  That  is,  for  i  <  if  e  Eq.  (11)  applies.  For  i  2:  if.t.  Eq. 
(13)  or  (14)  still  applies  for  the  static  temperature  at  the  reflection  line.  However,  Eq.  (12)  is 
modified  to  the  following  form: 


Pi.  1  v  i.  1  ~  2Pi.wvi,w  P i,2  V i , 2 

Pi.l  u.,l  =  fy1wui,w-Pi,2»l.2 

where  uiiW  =  -ViiW  sin  6\  and  vi  w  =  Vj(W  cos  0j  and  is  defined  by  Eq.  (15).  Substitution 
yields 

Pi  1  Vj  -1  =  2TTii  cosdi~  Pi,2Vi,2 

(34) 

Pi.l  ui.l  =  ~2mi  sin0i  “  Pi, 2  ui,2 

These  equations  allow  the  calculation  of  the  quantities  euvu  and  eu  Uj,i  but  the  value  of 
pi  !  is  needed  in  order  to  determine  Uj,j  and  Vj^.  This  is  accomplished  by  a  simple  iteration. 
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First  assume  that  Pi,i  =  pif2-  This  simulates  dp/dn  =  0  at  the  solid  wall.  Then  Tiit  is  known 
from  Eq.  (13)  or  (14).  Thus,  the  equation  of  state 


P.A  = 


. .  l 


T„ 


(35) 


gives  an  estimate  of  ^i.i  which,  with  the  use  of  Eq.  (34)  yields  uu,  vitI.  Now  Eq.  (17)  is  solved 
by  applying  Eq.  (18)  to  yield  an  improved  value  for  pu.  Then  Eq.  (35)  improves  eu.  This 
loop  is  repeated  as  necessary.  In  practice,  two  iterations  are  sufficient. 

4.2  MODIFICATION  TO  HYPERBOLIC  FINE  MESH  OPERATOR 


The  operator  which  solves  Eq.  (21)  governing  the  inviscid  contribution  in  the  fine  mesh  is 
called  L|J.  This  operator  was  originally  developed  based  on  the  assumption  that  V  <  <  a  in 
the  fine  mesh  where  a  is  the  sound  speed.  The  author  has  reservations  about  the  validity  of 
this  assumption  for  some  hypersonic  compression  corner  flows  with  separation  and  no  wall 
mass  transfer.  For  arbitrary  wall  mass  transfer,  this  assumption  is  clearly  invalid.  It  is 
therefore  necessary  to  reformulate  the  operator  Ljj.  Following  the  development  of 
MacCormack  (Ref.  2),  Eq.  (21)  is  expressed  as 


u. 


- 0 


(36) 


where  =  ( q ,  Uc,  Vc,  p)T,  Uc  =  u,  Vc  =  v  +  ijxu  and 

I  v..  Op  o 

Vp 


A  = 


P 

0 


0 

\ 0 


V 


yp 


\ 


■  / 


Uc  and  V,.  are  the  contravariant  velocity  components  in  the  computational  coordinates.  The 
compatibility  equations  for  this  hyperbolic  system  are  given  by 


<ip  pr 

<1  V 

< 

ilj'- 

ilp  pr 

c!  V 

l 

<U- 


1 1  nl 


<i.i- 


ulong  J—  =  i.  on.sLJ nt 


along  ,|  +■  =  ronsianl 


(37) 


(38) 
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where  c  =  aVl  +  and  J  +  and  J-  are  characteristic  coordinates  illustrated  in  Fig.  5.  These 
equations  are  identical  to  those  given  by  MacCormack  even  though  his  were  obtained  after 
assuming  that  Vc  <  <  a.  The  objective  for  the  operator  is  to  determine  the  time-averaged 
convection  and  pressure  fields  for  the  lime  step  of  interest  using  the  method  of 
characteristics.  The  conventional  MacCormack  scheme  is  then  applied  to  Eq.  (21)  using 


Figure  5.  Characteristic  coordinates. 

these  time-averaged  values  of  Vc  and  p  instead  of  the  local  values.  This  method  greatly 
enhances  the  stability  characteristics  of  the  operator.  This  operator  is  applied  at  a  givert_£ 
station  and  in  the  t-i\  plane.  Let  an  overbar  denote  a  time-averaging  operation.  Then  Pj,  VCj 
are  the  desired  quantities.  They  are  obtained  by  simply  integrating  their  gradients  to  obtain 


for  j  >  2  and 
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with  I'wnii  "5  I  '5  I’wl  and  V  =  Ap,  „ cost?,)  and  pw  is 

obtained  by  integrating  Eq.  (38)  along  the  I-  characteristic  to  the  wall.  The  quantities  dp/dri 
and  3Vc/3tj  are  obtained  from 

_  j_  r " +A<  i 

*,  =  ArJr,  rl,  '  (39) 

^  t  r-  Ar  ^vc  . 

fir,  =  a;  Jr  A?  ■  (40) 

where  dp/dtj  and  3Vc/3t|  are  expressed  in  terms  of  derivatives  with  respect  to  t  by  examining 
the  last  two  components  of  Eq.  (36)  given  by 


IV  ‘  YP  v ,  -  vJV,  =  0 

These  equations  may  be  expressed  as  QT  +  B  Q,  =  0  where 

0  =  (Vt  ,  p) 1  nntl 


Thus,  Q,  =  is  given  by  Q„  =  -B  'Qr  where 
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V2  -  c2 


By  direct  substitution,  Eqs.  (39)  and  (40)  become 


—  —r  +Ar  / 

dp  __  j_  f  °  l _ yp . 

dr,  At  Jto  \  V2  - 

7  _  j_  r»xA 7 _ 

dr,  ArJro  \  p(v2  ■ 


7  V 


(VC  - 


V2  -  c2 


V2  -  c2 


Pr  dr 


which  may  be  integrated  analytically  by  expressing  the  coefficients  of  VCr  and  pT  in  terms  of  a 
mean  value  over  the  interval  of  integration.  For  example,  in  the  general  case, 

f  A(x)  B(x)  dx  =  A(x)  /*  B(x)  dx  where  a  <  x  <  b 


provided,  of  course,  that  B(x)  is  of  one  sign  over  the  interval  of  integration.  The  result  is 


=  _L  y  p 

*7  Ar  v2  -  c2 
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'  +A t  V 

V  dr - i 

cr  v2 


V2  -  c2 


c  I 


,r  +At 


dr,  Ar  p  (y2  —  c2) 
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.ro  +  Ar 


pr  dr  =  p(rQ  -  Ar)  -  p(ro) 
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The  application  of  these  equations  at  point  j  gives 


where  p*  =  Pj(r0  +  At),  V*  =  vcj  (r0  +  At)  as  shown  in  Fig.  5  and  tj,  ti,  t3, 
r4,  have  been  taken  to  be  the  interval  midpoint.  The  unknowns  p*,  V*  are  determined  by 
solving  the  compatibility  equations  [Eqs.  (37)  and  (38)]  along  the  appropriate 
characteristics.  The  Finite-difference  representations  of  Eqs.  (37)  and  (38)]  are  given  by 


where  the  points  indicated  by  superscripts  P  and  M  are  illustrated  in  Fig.  5.  These  equations 
are  easily  solved  for  p*  and  V*  to  yield 


It  is  important  to  note  that  Eqs.  (41)  and  (42)  are  different  from  those  obtained  by 
MacCormack  (Ref.  2)  due  to  his  assumption  Vc  <  <  a.  In  particular,  the  terms  with 
coefficients  VC/(V2  -  c2)  were  neglected  by  MacCormack.  However,  these  terms  must  be 
included  in  the  operator  Ljj  because  of  the  new  wall  mass  transfer  boundary  condition. 
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The  next  section  provides  the  results  of  various  verification  phases  of  the  work  presented 
herein.  These  results  illuminate  many  of  the  same  characteristics  and,  in  some  cases, 
deficiencies  that  exist  for  the  original  two-dimensional  version  of  the  computer  code. 

5.0  RESULTS  AND  DISCUSSION 

All  of  the  numerical  results  presented  in  this  section  are  for  a  hollow  cylinder- flare 
configuration  in  a  Mach  2.98  freestream.  The  free-stream  static  temperature  is  202.85°R  and 
the  unit  Reynolds  number  is  10.97  million  per  foot. The  assumed  isothermal  wall 
temperature  is  taken  to  be  520°R.  The  configuration  is  28  in.  from  the  cylinder  leading  edge 
to  the  cylinder-flare  junction  point;  its  outer  radius  is  5.97  in.,  and  the  flare  angle  is  25  deg. 
The  flare  length  is  about  5  in.  This  configuration  was  tested  experimentally  by  Roshko  and 
Thomke  (Ref.  14).  The  numerical  results  of  the  present  study  are  compared  to  these 
experimental  results  wherever  possible. 

5.1  CODE  VERIFICATION 

An  initial  flow  solution  is  required  to  start  a  numerical  calculation.  Since  a  steady-state 
solution  to  the  Navier-Stokes  equations  is  sought,  this  initial  solution  is  arbitrary  within  the 
limitation  that  it  must  be  sufficiently  smooth  not  to  cause  failure  of  the  numerical  algorithm. 
The. computer  code  is  designed  to  provide  various  options  for  performing  this  initialization. 
One  of  these  options  is  used  when  a  solution  is  desired  over  an  entire  cylinder-flare 
configuration.  That  is,  the  computational  domain  extends  over  the  entire  length  of  the 
configuration  (see  Fig.  6),  The  leading  edge  singularity  is  captured  as  part  of  the  numerical 

—  —  —  —Computational  Domain  Over  Entire  Configuration 

- Computational  Domain  with  Boundary-Layer  Code 

Inflow  Solution 


solution.  This  option  is  used  when  the  interaction  region  at  the  cylinder- flare  corner  extends 
upstream  more  than  20  percent  of  the  distance  to  the  leading  edge.  In  this  case  the 
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interaction  region  is  adequately  described  by  several  grid  points.  If  this  interaction  extent  is 
much  less  than  20  percent  of  the  distance  to  the  leading  edge,  then  too  few  grid  points  are 
available  to  describe  it  and  the  flow  field  over  the  entire  cylinder.  In  this  case,  a  second 
option  is  used  for  initialization.  That  is,  a  solution  is  obtained  at  some  x-location  such  as  x*, 
well  ahead  of  the  interaction  region  using  a  boundary-layer  code  attributable  to  Patankar 
and  Spalding  (Ref.  22).  This  solution  is  then  used  to  define  an  inflow  boundary  condition 
for  the  new  computational  domain  which  has  this  station,  x*,  as  its  left-most  boundary.  This 
option  was  used  in  the  present  study  since  the  experimental  data  of  Roshko  and  Thomke 
(Ref.  14)  indicate  an  interaction  length,  f,  of  about  1  in.  Thus  the  interaction  extent  is  about 
3.6  percent  of  the  distance  from  the  cylinder-flare  corner  to  the  cylinder  leading  edge.  If  the 
first  initialization  option  were  used,  approximately  200  grid  points  would  be  required 
between  the  cylinder  leading  edge  and  the  end  of  the  flare  to  provide  only  six  points  within 
the  interaction  region  (assuming  equal  grid  point  spacing  in  the  x-direction).  This  is  clearly 
an  enormous  waste  of  computational  effort  when  the  interest  is  in  the  interaction  and  its 
effect  on  the  flare  solution. 


Figure  7  illustrates  the  actual  finite-difference  grid  resulting  from  the  use  of  the  second 
initialization  option.  The  abscissa,  x't  is  measured  from  the  cylinder-flare  corner,  and  the 
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The  interaction  length  is  defined  to  be  the  distance  from  the  beginning  of  the  pressure  rise  to  the 
corner  (see  Fig.  8b). 
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ordinate,  y,  is  measured  from  the  centerline.  Note  the  heavy  concentration  of  grid  lines  near 
the  wall  boundary  to  resolve  the  viscous-dominated  flow  features. 

As  part  of  the  verification  of  the  axisymmetric  capability  of  the  computer  code,  two 
numerical  solutions  were  obtained  at  zero  flare  angle.  These  two  solutions  were  obtained 
with  the  “current”  turbulence  model  of  Cebeci,  Smith,  and  Mosinskis  (Ref.  20)  and  the 
“alternate”  turbulence  model  of  Baldwin  and  Lomax  (Ref.  21)  (see  Section  2.3  for  a 
description  of  these  models).  The  axial  velocity  profiles  from  these  two  solutions  at  station 
x'  =  0  (obtained  with  x*  s  0.8L)  are  compared  in  Fig.  8a  to  the  profile  obtained  from  the 
Patankar-Spalding  boundary-  layer  code  at  the  same  streamwise  station.  The  solution 
obtained  with  the  “alternate”  turbulence  model  shows  better  agreement  with  the  Patankar- 
Spalding  profile  than  the  solution  obtained  with  the  “current”  turbulence  model.  For  this 
reason  and  because  the  boundary-layer  code  has  been  verified  against  experimental  data  by 
Patankar  and  Spalding  (Ref.  22),  all  subsequent  calculations  were  made  using  the 
“alternate”  turbulence  model. 

Additional  verification  of  the  axisymmetric  capability  of  the  computer  code  is  illustrated 
in  Fig.  8b.  A  comparison  of  surface  pressure  distributions  on  the  actual  hollow  cylinder- 
flare  configuration  is  shown  here  for  the  numerical  solution  and  the  experimental  data  of 
Roshko  and  Thomke  (Ref.  14).  The  computed  pressure  rise  and  asymptotic  value  on  the 
flare  section  agree  well  with  the  experimental  data.  The  interaction  extent  for  the  numerical 
solution  is  about  4um  =  0.5  in.  while  the  experimental  results  suggest  an  extent  of  about  4xp 
=  1  in.  Although  this  discrepancy  is  significant  when  the  details  of  separation  are  of 
interest,  it  has  nothing  to  do  with  the  axisymmetric  modification  to  the  computer  code.  This 
comment  is  strongly  supported  by  strikingly  similar  comparisons  which  are  documented,  for 
example,  in  Refs.  8,  9,  and  10  for  a  number  of  two-dimensional  problems.  These  references 
also  indicate  that  improvements  in  the  numerical  solution  result  from  the  application  of 
special  ad  hoc  features  to  the  turbulence  models  such  as  the  so-called  relaxation  approach.  It 
is  beyond  the  scope  of  this  effort,  however,  to  experiment  with  such  features. 

It  must  be  pointed  out  that  the  largest  uncertainty  in  the  state-of-the-art  numerical 
solutions  to  high  Reynolds  number  turbulent  flows  is  attributable  to  the  failure  of  the 
turbulence  modeling  effort  to  accurately  reflect  the  physical  turbulence  characteristics  of  the 
flow.  A  great  deal  of  research  is  currently  being  devoted  to  this  turbulence  modeling 
problem  worldwide  because  of  the  important  role  of  turbulence  in  modern  aerodynamics 
design  work . 
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a.  Velocity  profiles  at  x'=  0  for  zero  flare  angle  condition 
Figure  8.  Test  of  computer  code  on  a  simple  axisymmetric  turbulent 
boundary  layer. 


b.  Cylinder-flare  surface  pressure  distribution 
Figure  8.  Concluded. 
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5.2  FLOW-FIELD  DETAILS 

This  section  illustrates  some  details  of  the  shock-induced  separated  flow  field  with  the 
use  of  Figs.  9a  through  d.  Figure  9a  shows  contours  of  stream  function  which  clearly 
illustrate  the  recirculation  zone  and  the  flow  lift-off  occurring  near  the  separation  point.  The 
wiggles  appearing  in  these  contours  are  caused  by  numerical  “noise”  in  the  solution  in  the 
vicinity  of  the  compression  corner.  These  wiggles  appear  in  spite  of  the  use  of  a  product 
fourth  order  artificial  viscosity  term  by  Hung  and  MacCormack  (Ref.  1).  It  is  important  to 
note  that  the  recirculation  bubble  serves  as  a  sharp  turn  smoothing  device  which  converts  an 
abrupt  compression  corner  into  a  gradual  compression  followed  by  a  gradual  expansion 
followed  by  another  less  gradual  compression.  The  streamlines  over  the  separation  bubble 
illustrate  this  behavior. 

An  illustration  of  the  velocity  profiles  showing  the  onset  of  reversed  flow  (separation) 
and  reattachment  is  given  in  Fig.  9b.  Note  the  full  profile  at  x'  =  -0.45  in.  In  contrast  to 
this  profile,  the  profile  at  x'  =  +0.6  in.  behind  the  separation  zone  is  not  full.  In  fact,  the 
flow  does  not  develop  into  a  full  profile  behind  the  separation  zone  until  it  reaches  about 
half  way  back  on  the  flare  (i.e.,  x'  as  2.0  in.). 

The  shock  siructure  is  illustrated  by  static  pressure  contours  (isobars)  in  Fig.  9c.  The 
shock,  which  is  a  physical  discontinuity  for  all  practical  purposes,  is  shown  to  be  smeared 
out  in  this  figure.  This  is  attributable  to  the  application  of  the  so-called  shock-capturing 
approach  in  which  no  special  provision  is  made  in  the  computer  code  for  the  existence  of 
such  discontinuities.  If  they  exist,  they  are  “captured”  by  the  numerical  scheme.  A  key 
feature  of  this  figure  is  the  manner  in  which  the  contour  lines  all  approach  the  wall  in  a  near- 
normal  direction.  This  implies  that  dp/dn  =  0  at  the  wall  and  in  a  region  near  the  wall.  This 
is  consistent  with  boundary-layer  theory  which  is  perhaps  fortuitous  since  it  contrasts  with 
the  experimental  results  of  Settles,  Vas,  and  Bogdonoff  (Ref.  7)  who  show  contrary  results 
near  the  reattachment  region  for  a  two-dimensional  case. 

Figure  9d  is  a  simple  contour  of  the  flow  sonic  line.  This  figure  illustrates  how  abruptly 
the  subsonic  zone  grows  in  height  at  the  separation  point,  yet  its  retreat  back  near  the  wall  is 
asymptotic  and  requires  over  half  the  flare  length  to  reach  a  relative  height  equivalent  to  its 
height  just  ahead  of  separation.  This  is  consistent  with  the  fact  that  full  velocity  profiles  in 
the  boundary  layer  develop  slowly  after  reattachment  of  the  flow. 
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d.  Sonic  tine  contour 
Figure  9.  Concluded. 
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5.3  EFFECTS  OF  WALL  MASS  TRANSFER 

The  wall  mass  transfer  option  required  extensive  modification  to  the  computer  code  as 
mentioned  in  Section  4.0.  This  modified  code  was  first  verified  by  solving  the  hollow 
cylinder-flare  configuration  with  zero  wall  mass  transfer.  These  results  were  expected  to 
compare  favorably  with  the  numerical  solution  shown  in  Fig.  8b.  Figure  10  illustrates  that 
this  is  indeed  true.  In  fact,  the  code  with  the  wall  mass  transfer  modifications  produced  a 
solution  which  is  slightly  closer  to  the  experimental  data  of  Roshko  and  Thomke  than  the 
solution  from  the  code  without  these  modifications.  The  slight  difference  between  the  two 
numerical  solutions  is  caused  by  the  modifications  made  to  the  hyperbolic  fine  mesh 
operator  and  to  the  wall  boundary  conditions  described  in  Section  4.0. 

The  next  test  was  to  introduce  suction  into  the  cylinder-flare  corner  and  qualitatively 
analyze  the  results.  The  wall-normal  mass  flux  for  this  test  was  taken  to  be  -  0.01  gcue  where 
6eUf  is  the  axial  mass  flux  at  the  edge  of  an  undisturbed  boundary  layer  28  in.  long  (the 
length  of  the  cylinder  portion),  and  the  negative  sign  indicates  suction  rather  than  blowing. 
This  normal  flux  was  applied  at  one  point  on  each  side  of  the  cylinder-flare  corner  (there  is 
no  point  exactly  at  the  corner),  giving  an  axial  suction  extent  of  0.13  in.  on  each  side  of  the 
corner.  A  reduction  in  the  size  of  the  separation  region  was  expected  from  this  suction. 
Figure  11  illustrates  that  about  a  30  percent  reduction  did  result,  and  the  intermediate 
pressure  plateau  which  appeared  between  -0.25  in.  <  x'  <  0.25  in.  for  the  zero  wall  mass 
transfer  case  vanished.  The  solution  away  from  the  corner  is  unaffected  by  the  suction.  As 
mentioned,  this  test  is  qualitative  in  nature.  Further  testing  of  this  option  is  needed  for  a 
laminar  case  for  which  experimental  data  are  available  and  where  the  uncertainty 
attributable  to  turbulence  modeling  is  eliminated. 

6.0  CONCLUDING  REMARKS 

The  computer  code  documented  in  Refs.  1  and  2  by  Hung  and  MacCormack  was 
extended  to  apply  to  axisymmetric  flows  over  cylinder-flare  shapes.  Two  new  operators 
developed  for  this  extension  were  added  to  the  existing  time-split  sequence  of  operators.  The 
operator  L^,  which  involves  radial  derivative  term  contributions  to  the  complete 
axisymmetric  source  term  group,  was  found  to  have  an  upper  stability  bound  on  the  step 
size,  At,  which  is  governed  by  the  CFL  condition.  This  bound  is  proportional  to  the  Prandtl 
number  and  inversely  proportional  to  the  ratio  of  thermal  capacities  of  the  fluid,  7,  times  the 
maximum  value  over  ail  grid  points  of  the  kinematic  viscosity  scaled  by  the  radial  coordinate 
(see  Section  3.3).  It  was  found  in  practice  that  this  stability  bound  was  several  orders  of 
magnitude  less  restrictive  for  the  cylinder-flare  problem  examined  in  this  study  than  bounds 
already  present  attributable  to  the  other  operators  in  the  sequence;  however,  this  may  not 
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Figure  10.  Verification  of  new  wall  mass  transfer  boundary 
condition  and  associated  operator  for  the  case 
of  zero  wall  mass  transfer. 
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Figure  11.  Effect  of  suction  at  the  cylinder-flare  junction. 
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always  be  the  case.  Two  algebraic  eddy  viscosity  turbulence  models  were  described  and 
tested.  The  model  of  Cebeci,  Smith,  and  Mosinskis  (Ref.  20)  was  discarded  in  favor  of  the 
Baldwin-Lomax  model  (Ref.  21)  based  on  comparisons  made  with  a  boundary-layer  solution 
from  the  code  of  Patankar  and  Spalding  (Ref.  22).  Additionally  the  computer  code  was 
modified  to  allow  for  the  option  of  wall  mass  transfer.  This  required  development  of  a  new 
hyperbolic  Fine  mesh  operator,  Ljj,  because  of  the  assumption  in  the  development  of  the 
original  operator  that  V  <  <  c  in  the  fine  mesh.  The  final  computer  code  was  verified  by 
comparisons  of  boundary-layer  profiles  with  those  from  the  Patankar-Spalding  code  and  by 
comparisons  of  surface  pressure  distributions  with  the  experimental  results  of  Roshko  and 
Thomke  (Ref.  14)  for  a  particular  hollow  cylinder-flare  configuration.  The  computed 
surface  pressure  distributions  agree  poorly  in  the  separation  region  because  of  the 
inadequacies  of  the  turbulence  model,  which  is  entirely  consistent  with  computed  results  for 
two-dimensional  flows  over  compression  corners  (e.g.,  see  Refs.  8,  9,  and  10).  Flow-field 
details  are  presented  illustrating  shock  structure,  separation,  and  reattachment  details  as 
well  as  quick  buildup  and  gradual  decay  of  a  subsonic  pocket  near  the  cylinder-flare  corner. 
Qualitative  effects  of  introducing  suction  into  the  cylinder-flare  corner  are  presented  and 
discussed. 

The  final  computer  code  is  considered  applicable  in  general  to  flows  over  two- 
dimensional  compression  corners  and/or  cylinder- flare  configurations  at  zero  incidence.  In 
particular,  application  should  be  limited  to  laminar  flows  with  and  without  separation,  and 
turbulent  flows  for  which  an  eddy  viscosity  turbulence  model  adequately  describes  the  flow 
physics.  These  problems  may  have  suction  or  blowing  at  the  wall. 

Future  effort  is  needed  in  several  areas:  first,  the  turbulence  modeling  problem  is  by  far 
the  most  restrictive  and  should  therefore  be  given  the  most  attention;  second,  the  code  is 
designed  for  two  specific  geometric  shapes  —  wall-ramp  or  cylinder- flare  configurations. 
This  constraint  should  be  relaxed  so  the  flow  is  computable  over  virtually  any  two- 
dimensional  or  axisymmetric  shape;  finally,  the  effects  of  a  real,  rather  than  a  perfect,  gas 
should  be  included  to  extend  the  code’s  ability  to  provide  accurate  answers  to  real 
engineering  problems. 
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NOMENCLATURE 

A  Jacobian  matrix  used  in  stability  analysis 

A+  van  Driest ’s  constant 

a  Sound  speed 

B  Coefficient  matrix  used  in  development  of  new  Ljj 

b  General  constant  in  Klebanoff-type  intermittency  factor 

C  Constant  in  outer-layer  turbulence  model 

Cop  Constant  in  Baldwin-Lomax  turbulence  model 

Ckleb  Specific  Baldwin-Lomax  constant  for  Klebanoff-type  intermittency  factor 

C\yjc  Baldwin-Lomax  turbulence  model  wake  constant 

c  Geometry-scaled  sound  speed 

cp  Specific  heat  at  constant  pressure 

cv  Specific  heat  at  constant  volume 

D  van  Driest  damping  factor 

E  Flux  vector  used  in  discussion  of  time-splitting 

e,  e  Flux  vectors 

£  Specific  internal  energy 

F  Flux  vector  used  in  discussion  of  time-splitting 

f, f  Flux  vectors 

/  Denotes  functional  dependence 

G  Denotes  a  measure  of  local  mass-weighted  velocity  gradient  for  turbulence 

models 

g, g  Flux  vectors 

H  Total  mesh  height 
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Hf  Fine  mesh  height 

h,h  Source  term  vectors 

A 

i  Cartesian  unit  vector  in  x-direction 

A  A  A 

if,  in,  if  Curvilinear  unit  vectors  in  £,  ij,  f  directions,  respectively 

J  +  ,  J“  Characteristic  coordinates  in  the  fine  mesh 

a 

j  Cartesian  unit  vector  in  y-direction 

K  Constant  in  Cebeci-Smith-Mosinskis  turbulence  model 

k  Conductivity  coefficient 

k  Cartesian  unit  vector  in  z-direction  and  constant  in  turbulence  model  mixing- 

length  expression 

Li  Time-split  finite-difference  operator  associated  with  ^-direction 

L7  Time-split  finite-difference  operator  associated  with  ^-direction 

Time-split  finite-difference  operator  for  hyperbolic  terms  in  redirection 

Time-split  finite-difference  operator  for  parabolic  terms  in  ^-direction 

LhA  Time-split  finite-difference  operator  for  entire  source  term  with  conservative 

dependent  variables 

L*  Time-split  finite-difference  operator  for  entire  source  term  with  primitive 

dependent  variables 

L|  Portion  of  associated  with  ^-direction 

iJ  Portion  of  L*  associated  with  n-direction 

f  Mixing  length  or  separation  interaction  length 

M  Nonlinear  algebraic  transformation  matrix  relating  conservative  variables  to 

primitive  variables 

l 

m  Mass  flux 
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Pr  Prandtl  number 

p  Static  pressure 

Q  Klebanoff-type  intermittency  factor  and  dependent  variable  vector  used  in 

development  of  Lj| 

q  Velocity  vector 

s  Source  term  vector 

T  Static  temperature 

t  Time 

Uc  •  Contravariant  velocity  component  normal  to  £  =  constant  lines 
Udif  Parameter  used  in  Baldwin-Lomax  turbulence  model 

u  Axial  velocity  component 

V  Wall-normal  velocity  component 

Vc  Contravariant  velocity  component  normal  to  r)  =  constant  lines 

v  Velocity  component  in  y-direction 

w  Velocity  component  in  z-direction 

x  Axial  direction 

x'  Shifted  axial  coordinate 

x*  Location  of  left-most  grid  boundary 

y  Direction  normal  to  axis 

y  Reference  distance  normal  to  the  wall  for  turbulence  model 

y+  Inner-law  variable  used  in  turbulence  model 

z  Cartesian  coordinate 

7  Specific  heat  ratio 


51 


AEDC-TR-80-24 


Ar  Temporal  step  size 

A£,  At)  Spatial  computational  mesh  increments 

&  Boundary-layer  thickness 

4*  Displacement  thickness 

e  Total  energy  per  unit  volume 

t >V  Curvilinear  coordinates  in  an  axis-normal  plane 

8  Ramp/flare  angle 

*  Bulk  viscosity  coefficient 

X  Second  coefficient  of  viscosity 

p  First  coefficient  of  viscosity 

v  Kinematic  viscosity 

v '  Second  coefficient  of  viscosity  divided  by  density 

£  Computational  axial  coordinate 

q  Density 

o  Eigenvalue 

r  Computational  temporal  coordinate 

Tr  Reference  shear  stress 

$  Solution  vector  used  in  discussion  of  time-splitting 

<t>,  0  Solution  vectors 

0  Primitive  variable  vector 

0'  Reduced  primitive  variable  vector 

w  Vorticity 

V  Gradient  operator 
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SUBSCRIPTS 

B  Indicates  body  or  wall  value 

c  Indicates  corner  value 

e  Denotes  edge  value 

i  Denotes  inviscid  or  represents  index  in  x-direction 

j  Represents  index  in  y-direction 

f  Denotes  laminar  quantity 

f.e.  Denotes  leading  edge 

r  Denotes  reference  value 

t  Denotes  turbulent  quantity;  also  indicates  partial  differentiation  with  respect  to 

the  independent  variable,  t 

v  Denotes  viscous 

w  Indicates  wall  value 

x  Indicates  partial  differentiation  with  respect  to  x 

y  Indicates  partial  differentiation  with  respect  to  y 

z  Indicates  partial  differentiation  with  respect  to  z 

f  Indicates  partial  differentiation  with  respect  to  f 

■q  Indicates  partial  differentiation  with  respect  to  rj 

£  Indicates  partial  differentiation  with  respect  to  £ 

t  Indicates  partial  differentiation  with  respect  to  r 

SUPERSCRIPTS 

a,  a  Fractional  step  indicators 

b,  b  Fractional  step  indicators 

c  Denotes  corrected  value 
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h  Denotes  hyperbolic 

M  Denotes  minus  for  the  method  of  characteristics 

m  General  fractional  step  index 

n,  n  +  1 ,  n  +  1  Fractional  step  indicators 
P  Denotes  plus  for  the  method  of  characteristics 

p  Denotes  parabolic  or  predicted  value 

T  Denotes  transpose  operation 

y,  £  Indicate  source  term  groupings 

*  Denotes  solution  point  for  the  method  of  characteristics 

OTHER 

(  )  Indicates  association  with  Cartesian  reference  frame  or  denotes  a  time- 

averaged  quantity 

(')  Indicates  association  with  cylindrical  coordinates  (except  for  unit  vector 

and  other  specifically  mentioned  applications) 

(  )•(  )  Indicates  dot  product 

(  )x(  )  Indicates  cross  product 
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